Simplistic Coulomb forces in molecular dynamics: Comparing the Wolf and 

shifted-force approximations 
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This paper compares the Wolf method to the shifted forces (SF) method for efficient computer 
simulation of isotropic systems interacting via Coulomb forces, taking results from the Ewald sum- 
mation method as representing the true behavior. We find that for the Hansen-McDonald molten 
salt model the SF approximation overall reproduces the structural and dynamical properties as ac- 
curately as does the Wolf method. It is shown that the optimal Wolf damping parameter depends on 
the property in focus, and that neither the potential energy nor the radial distribution function are 
useful measures for the convergence of the Wolf method to the Ewald summation method. The SF 
approximation is also tested for the SPC/Fw model of liquid water at room temperature, showing 
good agreement with both the Wolf and the particle mesh Ewald methods; this confirms previous 
findings [Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006)]. Beside its conceptual simphcity 
the SF approximation implies a speed-up of a factor 2 to 3 compared to the Wolf method (which is 
in turn much faster than the Ewald method). 
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I. INTRODUCTION 



In molecular dynamics simulations the force evaluation 
consumes by far the most computational resources. For 
relatively short-ranged interactions like van der Waals in- 
teractioE^i it is common to introduce a cutoff radius Vc 
such that if the distance between a particle pair exceeds 
Tc, the particles do not interact^. This truncation allows 
for different optimization methods like inclusion of cell 
and neighbor lists^^^, which increase computational per- 
formance considerably. Traditionally, the pair potential 
is simply truncated and shifted such that it is zero at 
r^^—. This does not affect the force acting between par- 
ticles at distances below Tc, and if Vc is sufficiently large, 
the fluid properties are virtually unaffected by this ap- 
proximation. In fact, it has been showni*^ that keeping 
merely the short-ranged, purely repulsive part of the van 
der Waals interaction can account for the fluid structure 
even near the critical point where correlations are long 
ranged. The truncated and shifted potential approxi- 
mation ensures continuity of the potential energy, but 
introduces a discontinuity in the force at Tc, leading to 
energy drift for long simulation times^. To overcome this 
one can instead apply a truncated and shifted force (SF) 
approximation^, which has superior numerical stability®. 
Beside the numerical stability, it was recently shown by 
Toxvaerd and Dyre^ that for highly dense fluids the SF 
method allows for very small cut-off radius of rc ~ 1.5a- 
(where a is the atomic diameter) and corresponds ap- 
proximately to the first local minimum in the radial dis- 
tribution function. Applying such low cutoff to the trun- 
cated and shifted potential will lead to wrong physics 
and large energy drift'^ . The SF method therefore de- 
creases the number of interactions significantly and thus 
the computational time. The potential corresponding to 
the SF interaction does however not match the original 
potential for r < Vc from which the SF interaction was 
derived. Therefore, the thermodynamical properties can- 



not be compared directly, but can be derived from per- 
turbation theorySiii. 

For long-ranged interactions, like the Coulomb inter- 
action, one cannot simply introduce a standard cut and 
shifted potential. For example, simply truncating and 
shifting the Coulomb potential produces spurious fluid 
structure and wrong dynamics^. Numerous attempts 
have been made to overcome this problem. For example 
it has been suggested to use smoothing functions, but 
this leads in general to poor results, see Refs. [olliol Wolf 
et alM cleverly showed that using a simple truncated 
and shifted Coulomb potential corresponds in practice to 
summing over interactions in a non-neutral sphere. To 
compensate for this these authors introduced a neutraliz- 
ing term into the Coulomb potential; they further showed 
that faster convergence to the true energy is achieved by 
applying a damping factor a. The Wolf method is com- 
putationally much faster than the classical Ewald sum- 
mation technique and is today widely used within the 
scientific simulation community. The choice of the damp- 
ing factor, a, is, like the Ewald damping parameteii^ii^, 
somewhat arbitrary, and the optimal value must be found 
by comparison with either experimental data or results 
from, e.g., the Ewald metho d^"'^'^ . If the Wolf damping 
parameter a is zero, the Wolf method reduces to the SF 
approximationii, see also Denesyuk and Week&i^ for a 
discussion. We note that an SF method for the Coulomb 
interactions was used as a clever trick in the biochemi- 
cal simulation communityi^iii before the work by Wolf 
et at. 

In this paper we apply the Wolf method in molecular 
dynamics simulations of a simple model of a molten salt 
and liquid water. In order to find the optimal value of 
the Wolf damping parameter a we compare the simulated 
thermodynamical, dynamical, and structural properties 
with previously published results^^ based on the Ewald 
method. We show that the optimal value of a depends 
on the property one wishes to calculate and the cutoff 
distance used. This sets the stage for documenting the 
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main conclusion of this paper: for the systems studied 
here the SF approximation works as well as the Wolf 
method, confirming similar findings of Fennell and Gezel- 
ter— . Besides being conceptual simpler than the Wolf 
method, the SF method allows for more than a doubling 
of the computational speed. 



II. THE WOLF APPROXIMATION TO THE 
COULOMB POTENTIAL 

If r is the distance between two particles, the force 
acting on one particle from the other is F(r) = /(r)r/r, 
where / is for simplicity denoted the "force" and is mi- 
nus the derivative of the corresponding potential func- 
tion with respect to r. For the Wolf method"'^^ the force 
is given by 



Jw[r;a,rc) 



erfc(ar) erfc(c 



2J2\ „,2^2^ 



2a / exp(— a^r^) exp(— a^r^) 



(1) 



for r < Tc and where erfc(a;) — 1 — erf(x) is the comple- 
mentary error function. Here Zi and Zj are the charges 
of the two particles in question, Tc is the cutoff (i.e., 
fw = for r > Tc), and a is the Wolf damping pa- 
rameter. In the paper by Wolf et al. it is implicitly 
understood that arc > 1 such that the cutoff only takes 
effect beyond range of damping. The damping parameter 
was introduced in order to ensure faster convergence to 
the limiting Madelung energy^i. Unfortunately, there is 
no theoretical prediction for the optimal value of a, which 
must be found by comparison with other well-established 
methods like the Ewald summation method^'^'^^'^'^. Wolf 
et ali^ and Demontis et al^ have shown via molecular 
dynamics simulations that the Wolf method reproduces 
the results obtained by the Ewald summation method for 
Tc > 5(iij, where dij is the distance between oppositely 
charged particles in the first coordinate shell. Demontis 
et ali^ also suggested that the optimal damping param- 
eter is given by a = 2/rc for sufficiently large systems. 

From Eq. ([TJ it follows that for a — ^ oo one has fw 
0, and that for a — > the force reduces to 



fspir; rc) = z^zj (l/r^ - l/r^) for r < 



(2) 



This is the truncated and shifted force (SF) cutoff^'^. 

In Fig. [IJa) we plot the difference between the Wolf 
force, fw, and the corresponding Coulomb force, fc = 
ZiZj/r^ , for different damping parameters. Clearly the 
damping parameter has a non-trivial effect on the force. 
For a = the difference is small compared to large val- 
ues of a, suggesting that the SF cutoff, Eq. ([2]), gives a 
good approximation to the Coulomb interaction. From 
Fig. [IJa) it is seen that an optimal value of a exists that 
minimizes the difference. One way to identify this opti- 
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FIG. 1: [Color online] (a): Difference between the Wolf force, 
/w, and the Coulomb force, fc, for a = 0.0, 0.25, 0.5, 0.75 and 
1.0. In all graphs the cutoff is given by = 4.18. (b): The 
measure of the difference between the true Coulomb force and 
the Wolf force, Ef, defined in Eq. ([3| plotted as a function 
of a for three different cutoffs. The inset shows the optimal 
value of a plotted as a function of rc. 



mal value is by minimizing the function 

Ef{a,rc) = 1 



fw{r;a,rc)dr 

° , (3) 

fc{r) dr 



which measures the total relative difference between 
fw{r) and fc{r) such that Ef > (since fw < fc for 
all r). In Fig. [ijb) Ef is plotted for three different 
cutoff distances. The optimal Wolf damping parameter 
converges to zero as rc increases, which reflects the sim- 
ple fact that fw fc for Tc — oo and a — 0. More 
interestingly, the quantity Ef does exhibit very little dif- 
ference between the optimal value of a and a = 0. The 
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inset in Fig. [Ub) shows that the optimal Wolf parameter 
determined by the minimum of Eq. ([3]) is given roughly 
by a w 3/(4rc). This simple analysis is consistent with 
the Tc dependence suggested by Demontis et al}^ based 
of molecular dynamics simulations (but they predict a 
smaller estimate of a by a factor of 3/8). 

The conclusion from Fig. [T] is that setting a = 0, 
i.e., adopting the SF approximation, gives results that 
are close to those obtained by carefully optimizing a. 
This and the recent work by Toxvaerd and Dyre'^ moti- 
vate the below reported molecular dynamics simulations, 
which compare the Wolf method to the SF cutoff for other 
quantities and realistic systems. As "truth" we take the 
well-established, but computationally expensive, Ewald 
summation method. 
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III. RESULTS FOR THE HANSEN AND 
MCDONALD MOLTEN SALT MODEL 

A series of molecular dynamics simulations was per- 
formed of a model molten salt proposed by Hansen and 
McDonaldi^. Briefly, in this two-component model the 
ions are simple spherical particles that interact via a 
Coulomb potential and a van der Waals type potential 
given by the inverse power law <^(r) = ^ where 
n = 9, e defines the energy scale and a is the usual 
Lennard-Jones length scale parameter^. We refer the 
reader to the reference for the full details. In the sim- 
ulations we applied the Wolf method and varied the cut- 
off between 2.5 and 8.0 a. The simulation box used was 
twice the size of the cutoff whenever Tc > 4.18(7; for 
smaller cutoffs the box length was fixed to 8.36 a. The 
number density for all systems were p = 0.368(7^'^, thus, 
the number of ions varied from 216 to 1508. The re- 
sults presented below were found to be independent of 
system size. The temperature T is controlled using a 
Nose-Hoover thermostat^^-^^ with T = O.OlTTe/A:^. The 
results are compared to previously published data where 
the Ewald summation method was used^^, which repre- 
sent the "true" Coulomb interaction. 

First, in Fig. [5] (a) we compare the total potential en- 
ergy obtained from the Wolf method Uw for three differ- 
ent cutoff radii and varying damping parameters with the 
potential energy Ue from the Ewald summation method. 
We note that Uw is obtained directly from the Wolf po- 
tential functioi*i2iiii corresponding to the force given in 
Eq. dD). It is observed that Uw is within the statistical 
uncertainty equal to Ue for sufficiently small damping 
parameters, even for quite small cutoffs. This could lead 
to the conclusion that the Wolf method accounts cor- 
rectly for electrostatic interactions for small cutoff dis- 
tances. However, if one plots the radial distribution func- 
tion g, Fig. [51^b), we see that for Vc — 2.5a the structure 
differs significantly from the result obtained using the 
Ewald summation method. This is true for all values of 
the damping parameter. From Fig. [^b) we also notice 
that the SF approximation captures the structural prop- 




FIG. 2: [Color online] (a): Comparison of the potential energy 
of the Hansen-McDonald molten salt model for varying damp- 
ing parameter. Error bars represent the standard error of ten 
independent runs. Ue is found in Ref . [l^. (b): Radial distri- 
bution functions for unlike charged particles (lines) for a = 
(the SF approximation) and rc = 6a and for a = OAa^^ and 
Vc = 2.5a. The filled black squares are data points taken from 
Ref. M- 

erties correctly for rc = 6(t, which is the smallest cutoff 
distance meeting the Wolf et aL— and Demontis et ali^ 
criterion rc > 5dij. 

We study the radial distribution function dependence 
of Tc and a by defining the error parameter Eg via 

/ \9w{r) - gE{r)\dr 
E, = ^ j^^ , (4) 

Jo 

where gw is the radial distribution function for unlike 
charged particles of the Wolf method and gE the ra- 
dial distribution function produced by the Ewald sum- 
mation method. Similarly, the following error parameter 
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FIG. 3: [Color online] Error parameters as a function of cutoff 
for different damping parameter for the Ifansen and McDon- 
ald molten salt system. The inset shows the error parameters 
for rc = 8a as functions of a. 

Ed quantifies the difference in diffusion constant 

Ed = ^-1. (5) 

where Dw and De are the diffusion constants obtained 
from the Wolf and Ewald methods, respectively. Note 
that Eg > 0, whereas Ed can be negative. The "correct" 
radial distribution function, qe,, and diffusion constant, 
De, were taken from Ref. [l^. Figure |3] shows the two 
error parameters for different cutoff radii and damping. 
The damping parameter a = 0.7a~^ was chosen because 
Eg exhibits a minimum for this value for a large range of 
cutoffs. This is not the case for Ed, however, which fea- 
tures a minimum for lower values of the damping parame- 
ter, depending on the cutoff (as expected from Fig.[T](b)). 
This inconsistency is illustrated in the inset in which the 
error parameters are shown for rc = Su as functions of a. 
Obviously, any a < 0.6a^^ may be chosen to minimize 
Ed, whereas Eg features a minimum for a = 0.7a~^. We 
note that rc = 8a > 5dij and the cutoff radius fulfills the 
criterion defined by Wolf et al. and Demontis et ai. 

From Fig. [3] it is seen that Eg is relatively large 
for small cutoffs (as expected), but that it for non-zero 
damping parameters quickly decreases and reaches al- 
most zero for rc > 4.0cr~^. For the SF approximation 
one needs rc > 6.0(7 in order to obtain the same accuracy 
in the radial distribution function. For large cutoffs the 
SF approximation results in better diffusion constants 
than the Wolf method with a — 0.7a~^. We could, 
of course, have optimized a with respect to the diffu- 
sion constant (giving a = 0.3a~^ for a large range of 
cutoffs). This, however, would decrease the agreement 
for the radial distribution function. This fact is high- 
lighted in Table I, where the error parameters are listed 



a [cr"'] Ed Eg 

0.0 (SF) 0.04 ± 0.02 0.0f7 ± 0.002 

0.3 0.03 ± 0.01 0.019 ± 0.002 

0.7 0.12 ± 0.02 0.010 ± 0.001 



TABLE I: Error parameters, Ed and Eg, for different values 
of the damping parameter, a = 0.3a~^ and a — 0.7a^^ 
correspond to the optimized values with respect to diffusion 
and radial distribution function, respectively, a = O.Ocr"^ 
corresponds to the SF approximation. 



for values of a optimized, respectively, with respect to 
the diffusion constant and the radial distribution func- 
tion (rc = 8.0(t). For comparison we also give the error 
parameters for the SF approximation. Within the statis- 
tical uncertainty there is no difference between the Wolf 
method using a = 0.3cr~^ and the SF approximation. 

Up to this point we have only discussed the structural 
and diffusive properties in the long time limit. To com- 
pare the short-time dynamics of the two methods we plot 
the velocity autocorrelation function Cvv{t) and the in- 
termediate scattering function in Fig. S] From Figs. [2] 
and [3] it was concluded that for small rc {rc ~ 4. Oct) and 
large a (a ~ 0.7) both the potential energy and the ra- 
dial distribution function are in excellent agreement with 
the Ewald summation method, but in Fig. 2] we clearly 
observe the short-time dynamics is not correct for this set 
of parameter values. This shows that the cutoff must be 
sufficiently large for the Wolf method to correctly account 
for all the fluid properties - but at such large cutoff the 
SF approximation may be applied instead since it results 
in the same accuracy. 



IV. RESULTS FOR A WATER MODEL 

We also tested the SF approximation for liquid wa- 
ter at the state point (T, p) = (300 K, 998 kg m"^) using 
the flexible single point charge (SPC/Fw) water modeP^. 
In this model the chemical bond and the bending an- 
gle are allowed to vibrate around their zero-force values. 
The model is easy to implement and has been shown 
to predict many bulk properties better than for exam- 
ple the SPG, SPC/E and TIP3P models2ii22. In Fig. 
[5fa) we plot the oxygen-oxygen radial distribution func- 
tion goo for the Wolf and SF methods. For comparison, 
data from Ref. are shown (filled squares), where the 
Coulomb interactions were evaluated using the particle- 
mesh Ewald (PME) method^^. The radial distribution 
function is reproduced reasonably well by both meth- 
ods. The SF approximation captures the liquid structure 
at least as well as the Wolf method, except at the first 
peak which is slightly underestimated. The radial dis- 
tribution functions for both the SF and the Wolf meth- 
ods are independent of the cutoff for radii larger than 9 
A, the value used by Zahn et al}'^; this corresponds to 
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FIG. 4: [Color online] (a): Normalized velocity autocorrela- 
tion function for the Wolf and the SF methods. Only the 
short time data are shown. The error bars are comparable 
to the size of symbols, (b): Coherent intermediate scattering 
function for wave-length k = 7.18(t~^ and Tc = 7. Oct. The 
horizontal line is the interpolated value of the static structure 
factor S{k) = F{k,0) taken from Ref. [l^. The time t is given 
in standard reduced molecular dynamics units. 

Tc ~ bdij since the oxygen- hydrogen distance is around 
1.8 A. In Fig. [SJb) the center-of-mass velocity autocorre- 
lation function is plotted for two different cutoffs for both 
methods. This dynamic property is largely independent 
of method and cutoff, as is the case for the liquid struc- 
ture. The same conclusion was reached by Fennell and 
Gezelter— . For the SF approximation we obtain a diffu- 
sion constant of 2.4 xlO~^ m^ s~^ and a shear viscosity 
of 0.78 xlO"^ Pa s. This can be compared with the 
experimental values 2.3 xlO~^ m^ s^^ and 0.85 xlO^"^ 
Pa s. It is also worth mentioning that Zahn et al~^ used 
a = 0.06(7""'^ in their simulations of (rigid) SPC/E water, 
but found that the potential energy was in better agree- 
ment with the Ewald method for even lower damping 



FIG. 5: [Color online] (a): Oxygen-oxygen radial distribution 
function for the SPC/Fw water model using SF and Wolf 
methods. The squares represent data taken from Ref. [2ll. (b): 
Normalized center-of-mass velocity autocorrelation functions. 
The inset shows a zoom of the time interval 0.065 to 0.32 ps. 
In both (a) and (b) cr = 3.16 A. 



parameters. For a — O.OGu ^ we observe no difference to 
the SF approximation. 



V. CONCLUDING REMARKS 

In conclusion, for simple molten salts and liquid wa- 
ter the SF approximation reproduces various properties 
as well as the Wolf method. The Wolf method has one 
more parameter than SF, and consequently this method 
may be optimized to give slightly better agreement with 
the Ewald summation method. Such an optimization, 
however, must be carried separately out for each prop- 
erty in focus and for each different system. Beside its 
simplicity (and thus easy-to-code feature), we found that 
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the SF approximation leads to a simulation speed-up of 
2-3 compared to the Wolf method. Of course, the ac- 
tual speed-up depends on the specific problem and the 
use of optimization techniques, but the calculation of the 
four terms in Eq. ([ij involves complicated mathematical 
functions and is deemed to consume considerably more 
computational resources than the simple SF approxima- 
tion. We wish to stress here that the paper of Wolf et 
al. was the first to correctly analyze why the SF approx- 
imation for Coulomb forces is superior to the standard 
truncated and shifted potential interaction model. 

Fennell and Gezelter^^ carefully analyzed an impres- 
sive number of different systems including simple crys- 
tals showing a good agreement between the SF method 
and the Ewald technique. In their conclusion the au- 
thors suggested that the SF approximation can also be 
used for confined geometries, thereby overcoming the en- 
forced periodicity in the unmodified Ewald method. We 



agree that the Ewald method can be problematic (even 
for periodic systems^i), but, the SF approach (as well as 
the Wolf method) is an approximation that suppresses 
the intrinsic long-ranged nature of the Coulomb interac- 
tions leading to an artificially molecular orientationSi^ 
in confinements. For confined fluids alternative methods 

have recently been adviced, see Refs. [iliili^. 
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